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WAVELET  DOMAIN  IMAGE  INTERPOLATION  VIA  STATISTICAL  ESTIMATION 

Ying  Zhu  Stuart  C.  Schwartz  Michael  T.  Orchard 

Electrical  Engineering,  Princeton  University 
Princeton,  NJ  08544 


ABSTRACT 

We  propose  a  new  wavelet  domain  image  interpolation 
scheme  based  on  statistical  signal  estimation.  A  linear 
composite  MMSE  estimator  is  constructed  to  synthesize 
the  detailed  wavelet  coefficients  as  well  as  to  minimize 
the  mean  squared  error  for  high-resolution  signal 
recovery.  Based  on  a  discrete  time  edge  model,  we  use 
low-resolution  information  to  characterize  local  intensity 
changes  and  perform  resolution  enhancement 
accordingly.  A  linear  MMSE  estimator  follows  to 
minimize  the  estimation  error.  Local  image  statistics  are 
involved  in  determining  the  spatially  adaptive  optimal 
estimator.  With  knowledge  of  edge  behavior  and  local 
signal  statistics,  the  composite  estimation  is  able  to 
enhance  important  edges  and  to  maintain  the  intensity 
consistency  along  edges.  Strong  improvement  in  both  the 
visual  quality  and  the  PSNRs  of  the  interpolated  images 
has  been  achieved  by  the  proposed  estimation  scheme. 


1.  INTRODUCTION 

Image  interpolation  involves  the  problem  of  resolution 
enhancement,  i.e.  recovering  a  high-resolution  image 
from  its  smoothed  version.  In  low-resolution  images,  the 
lost  high-frequency  components  associated  with  strong 
intensity  changes  are  the  main  information  to  recover. 
Strong  intensity  discontinuity  conveys  substantial  visual 
information  because  it  often  happens  around  various 
edges  produced  by  object  boundaries.  Therefore,  how  to 
process  the  intensity  changes  around  strong  edges  is  the 
central  issue  for  the  interpolation  problem. 

Classic  bilinear  and  bicubic  interpolation  methods 
impose  a  continuity  constraint  over  the  entire  image  and 
tend  to  produce  oversmoothed  edges.  They  also  implicitly 
assume  that  a  low-resolution  image  consists  of  samples 
from  its  high-resolution  version,  which  is  not  true  when 
the  imaging  sensor  has  a  spatial  averaging  characteristic 
over  the  sampling  lattice.  Edge  adaptive  interpolation 
schemes  [1,2]  perform  the  interpolation  in  selective 
directions  to  avoid  smoothing  across  edges,  with  extra 
efforts  made  to  deal  with  edge  related  artifacts  due  to 

This  work  was  partly  supported  by  the  ARO  grant  DAAH0496 1-0227  and  the  NSF  grant 
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imprecise  edge  positions  and  smoothness  assumed  for  the 
high-resolution  image.  Wavelet  based  methods  [3,4]  take 
advantage  of  local  signal  smoothness  (Holder  regularity) 
observed  through  multiple  scales  [5]  and  actively  enhance 
the  high  frequency  contents.  However,  the  general 
mathematical  conclusion  involved  is  obtained  by 
analyzing  continuous  signals.  When  interpolating  discrete 
time  signals,  ambiguities  about  the  locations  and  the  signs 
of  the  extrapolated  highband  coefficients  would  arise. 

In  this  work,  we  propose  a  new  wavelet  domain 
interpolation  scheme  that  explicitly  fills  in  important 
detailed  coefficients  to  recover  the  high-resolution  image. 
In  contrast  to  other  wavelet  domain  interpolation 
methods,  we  formulate  interpolation  as  a  statistical  signal 
estimation  problem,  i.e.  we  want  to  estimate  the  high- 
resolution  image  using  the  statistical  information  of  the 
low-resolution  signals.  We  also  characterize  edge 
behavior  by  a  parameterized  discrete  time  signal  to 
accurately  locate  edges  from  low-resolution  samples.  A 
linear  composite  minimum-mean-squared-error  (MMSE) 
estimator  is  proposed  to  solve  the  estimation  problem. 
The  composite  estimation  involves  a  parametric  edge 
model  and  local  image  statistics.  First,  local  edge 
behavior  is  determined  from  the  low-resolution  samples 
and  used  to  synthesize  the  detailed  coefficients.  Then,  a 
linear  estimator  minimizes  the  estimation  error  using  local 
statistical  information  of  the  enhanced  signals.  Both 
analysis  and  experiments  reveal  that  the  knowledge  of 
edge  behavior  and  local  image  statistics  enable  the 
composite  estimator  to  enhance  the  cross-edge  sharpness 
and  maintain  the  intensity  consistency  along  edges,  which 
are  essential  for  high  image  quality. 

In  the  following  sections,  we  start  by  characterizing 
edge  behavior  under  lowpass  filtering,  and  then  derive  the 
optimal  MMSE  estimator  for  image  interpolation.  Finally 
we  show  the  experimental  results  which  demonstrate  the 
strength  of  the  statistical  approach. 

2.  WAVELET  TRANSFORM  AND  EDGE 
BEHAVIOR 

Figure  1  illustrates  a  ID  discrete  wavelet  transform  with 
analysis  filters  ({/;(«)}, ]g(n)})  and  synthesis  filters 

({/t(n)}, {#(«)}) .  Discrete  time  signal  /  is  decomposed 
and  down-sampled  to  produce  two  subband  signals,  the 
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lowpass  scaling  coefficients  f  and  highpass  detailed 


Figure  1.  ID  Wavelet  transform 


coefficients  fd.  The  synthesis  process  gives  exact 
reconstruction  of  f  Wavelet  transform  on  2D  images  can 
be  realized  by  separate  ID  transforms  on  image  rows  and 
columns. 

To  characterize  edge  behavior  under  lowpass 
filtering,  we  adopt  (1)  as  the  underlying  continuous  edge 
model,  which  assumes  the  ideal  step  function  for  real 
object  boundaries  and  a  Gaussian  point  spread  function 
for  the  acquisition  system  [6], 

e(t;c,A,<j)  =  c  +  ^-(l  +  erf(-jL-J)  (1) 

2  V  2(7 

Here  erf  O')  =  (2 /  4tz)\ exp(-x2  /  2  )dx  .  With  a  normalized 
o 

sampling  rate,  the  corresponding  discrete  time  edge 
model  is  parameterized  as 

e(rv,c,A,o',d)  =  c  +  —  (\+erf(n  ))  (2) 

2  V2(7 

where  (A,  a,  d)  indicate  the  edge  strength,  transition  speed 
and  sampling  bias,  and  c  indicates  the  base  intensity. 

Let  a  discrete  time  edge  signal  f{n;c,A,a,d)  be 
smoothed  with  lowpass  filter  { h{n) }  and  downsampled  by 
2,  the  low  resolution  signal  f  still  appears  in  the  shape  of 
a  slightly  modified  edge  as  we  expect.  However,  if  we 
also  allow  the  detailed  signal  fd  to  be  produced  through 
highpass  filter  {g(u)l  and  downsampling,  we  observe  the 
important  fact  that  for  sharp  edges  with  small  a,  besides 
the  apparent  dependence  on  A  and  a,  the  local  detailed 
coefficients  in  fd  also  change  dramatically  and  nonlinearly 
as  the  sampling  bias  d  changes.  Figure  2  (b)-(c)  illustrates 
the  effect  with  c  set  to  1  and  b  set  to  0  and  0.5 
respectively.  Two  facts  about  edge  behavior  are  used  in 
the  following  analysis.  First,  a  strong  intensity  change  in 
the  high-resolution  signal  causes  a  strong  intensity  change 
at  the  corresponding  location  in  the  low-resolution  signal. 
Second,  the  knowledge  of  edge  parameters  (A,  a,  d)  are 
sufficient  for  determining  the  local  detailed  coefficients. 

Now  we  consider  /  to  be  the  unknown  ID  discrete 
time  signal  taken  from  an  entire  row  or  column  of  a  high- 
resolution  image  and  assume /can  be  approximated  as  the 
composite  signal  of  a  number  of  local  edges  at  different 
spatial  locations 

/(»)  =  c  +  y^-(l  +  erf(”  d'))  (3) 

'  2  V2fj; 

Suppose  that  only  the  smoothed  signal  /  is  available.  To 
recover  the  significant  detailed  coefficients  of  fd,  our  basic 


idea  is  to  obtain  an  estimate  of  local  edge  parameters  (A,-, 


(a) 

./;, .  L 

V  . 

•  •  ■  •  •  •  •  • 

i  * 

»  .  »  .  (b) 

fs  *  *  *  *  *  fd  , 

•  •  1 1  \  ,*  •  •  *  • 

(c)  • 

Figure  2.  (a)  Discrete  time  edge  model;  (b)-(c) 

Scaling  coefficients  /  and  detailed  coefficients 
[,.  (b)  o=l,  b=0\  (c)  o=l.  b=0.5. 

(T,,  d)  from  the  low-resolution  signal  f  and  use  this 
information  to  fill  in  the  important  coefficients  of  fd. 

3.  MMSE  ESTIMATOR  FOR  IMAGE 

INTERPOLATION 

We  formulate  image  interpolation  as  an  estimation 
problem,  i.e.  given  the  low-resolution  image  samples 
{fs(x,y)},  we  want  to  estimate  the  high-resolution  samples 
\f(x,y) } .  To  simplify  the  analysis,  here  we  only  consider 
the  problem  of  resolution  enhancement  in  one  direction. 
To  achieve  complete  resolution  enhancement,  we  apply 
the  interpolation  scheme  twice.  First  we  interpolate  the 
image  horizontally,  and  then  we  apply  the  same  scheme 
vertically  to  the  horizontally  interpolated  image.  Let  f, 
denote  the  to-be-interpolated  high-resolution  image  row 
with  vertical  index  y.  fsy  denotes  the  y_th  row  in  the  low- 
resolution  image  and  fdy  denotes  the  y_th  row  of  the 
unknown  detailed  coefficients.  Figure  3  helps  to  illustrate 
the  notation.  In  particular,  we  consider  the  estimation  of 
fix),  the  high-resolution  image  pixel  indexed  by  (x.y). 
Under  the  minimum-mean-squared-error  (MMSE) 
criterion  [7],  the  optimal  estimator  for/y(x)  is  given  by  the 
following  conditional  mean 

fy  (X)  =  E[fy  (X)  I  fs.  y  +  k;k  =  0+1,  •  •  •  ±  M  ]  (4) 

In  practice,  strong  correlation  exits  only  in  local  image 
region,  so  M  can  be  set  small.  Direct  evaluation  of  (4) 
requires  the  knowledge  of  joint  distribution  of  fy  and 
{fs,y+k}-  This  is  difficult  to  obtain  when/,  is  unknown. 
Consequently,  we  propose  the  following  linear  composite 
MMSE  estimator  as  a  replacement, 

-  M 

fy M  =  L  ak  ■  E[fy (a-)  I  fs  v+k]  +  C  (5 ) 

k--M 
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where  ak  and  C  are  scalars.  Denote  fy  =E[fy  \  fs,  y  +  /t] 
with  fy  (n)  =  E[fy(n )  I  fs,  y  + t-] ,  we  rewrite  (5)  as 

/,(*)=  Z  akfyk(x)+C  (6) 

k=-M 

The  estimator  is  implemented  in  two  steps.  First,  f,  is 
estimated  as  fy  using  the  information  from  the  k_th 

neighboring  row  fs,y+k  in  the  low-resolution  image.  We 
call  this  step  the  estimate  by  row.  It  involves  the  use  of 
ID  edge  model.  Then  a  linear  MMSE  estimator  estimates 
fy(x)  by  minimizing  the  mean-square-error(MSE)  over  the 

M  ~k 

linear  class  Z  ak  ■  f  (x)  +  C  .  Second  order  statistics  of 

k=-M 

fk  are  involved  in  finding  optimal  values  for  { ak }  and  C. 

3.1.  Estimate  by  row 

The  synthesis  off.  is  given  by 

fy={f,,y^2)*h+(fdty\2)*g  (7) 

where  T  2  denotes  the  linear  operation  of  upsampling  by 
2,  so  we  have 

fy  =  (fs,y  t  2 )*h+  (. E[fd  y  I  fs  y+k]  t  2 )*g  (8) 

We  apply  the  ideas  discussed  in  section  2  and  solve  the 
estimate  fdy  =  E[fd  y  I  fs^y+k] .  Assume  the  k_th 

neighboring  row  fy+k  can  be  approximated  as  the 
composite  signal  of  a  number  of  edges 
c  +  Z;  A,(l  +  erf((n  -  df/  f2<j.))/2 .  Around  dj  there  is  an 

edge  segment  crossing  row  fy+k  that  causes  the  local 
intensity  change  on  fy+k  as  fy+k  (, n )  =  e(n;  c,. ,  A,. ,  er, ,  dt ) . 

Then  the  k_th  neighboring  row  f,y+k(n)  in  the  low- 
resolution  i mage  /’  has  a  local  intensity  change  around  d, 
/2  as 

Z,„  him)  ■  e(2n  -  m;  cj ,  A ,  er ,  dj )  (9) 

So  the  parameters  (c, ,  Al ,  er  ,  dt )  can  be  estimated  from 
fSjy+k  around  d,  /2  by  solving  the  least  square  (LS)  problem 
argmin  I  I  fty+k  (n)  -Zh(m )e(2n  -  m;c: ,  At , a,  ,dt)\2 

(c.A.d.cr)  neN(dj/2)  ’  m 

(10) 

Once  the  LS  estimates  (c. ,  ,  os ,  dt )  are  obtained,  we 

estimate  the  direction  of  the  edge  segment  crossing  row 
fy+k.  As  illustrated  in  figure  3,  we  use  the  edge  segment 
fs,y+k+j  (n)  =  Zm  h(m)  ■  e(2n  -  m;  c. ,  A,. ,  ai ,  di  -  j  /  tan(6> ))  (1 1) 

to  match  the  local  pixels  around  ( d/2,y+k )  in  f  and 
determine  the  edge  orientation  0,.  Then  the  estimates  of 
detailed  coefficients  around  d/2  in  fdy  is  given  by 

fly(.n)  =  Ztn,g(m)-e(2n-m;ci,Ai,ai,di+k/tan(0i))  (12) 

The  estimate-by-row  step  is  summarized  as  follows: 

1.  Check  the  k_th  neighboring  row  f,y+k  of  the  low- 
resolution  image.  If  a  strong  intensity  change  is  found 


/^i  dj+k/tan(8j) 

1  <5- fy-l 


fs.y-1  ( fd.y-1 V  y- 

fs.y(fd.y)  y  fy 

fs.y+1  (fd,y+l)* fy+] 
fs.y+k  (fd.y+k )  fy+k 

Figure  3.  Interpolation  lattice 


at  some  location,  perform  the  LS  estimation  in  (9)  to 
obtain  the  local  edge  parameters  (c. Use 
(10)  to  estimate  edge  orientation  0, 

2.  Use  (12)  to  synthesize  the  detailed  coefficients  in  fd 

around  d, 12  .  Then  use  (8)  to  obtain  the  estimate  fy  . 

3.  Apply  the  same  procedure  to  each  row  of  the  image. 
For  every  row  f,  we  obtain  (2M+1)  high-resolution 
estimates-by-  row,  fy  =  E[fy  I  fSjy+k]  (k=-M,...M). 

Based  on  the  family  of  parametric  edge  models,  the 
estimate-by-  row  implements  a  nonlinear  estimator  for  the 
high-resolution  signal  using  partial  information  from  the 
low-resolution  signal.  As  we  will  see,  resolution 
enhancement  around  strong  edges  is  achieved  by  the 
estimates. 


3.2.  Linear  MMSE 


The  linear  MMSE  (6)  finds  the  optimal  linear 
combination  of  the  estimates-by-row  that  minimizes  the 
estimation  MSE.  Only  first  and  second  order  statistics  of 
fy  are  involved  to  solve  the  optimal  {a*}  and  C  [7], 


M 


C  =  E[f  (x)]  -  £  ak  ■  E[fy(x)] 

k=-M  * 


(13) 


Cov(f  (x),fy(x))=  Z  a  :Cov(fy  (x),  fy  (x)) 


j=-M 


where  Cov(y )  denotes  covariance  function.  The  following 
facts  are  easily  verified: 


Elfy(x)]  =  E[fv(x)] 


COV(f  (X),  fy  (X))  =  COVijy  (X),fy  (X)) 


(14) 


Denote 


CUy) (;,  k)  =  CoviJ'y  (x),  /;  (x)) .  (14)  allows 


{ak}  and  C  to  be  solved  by  2M+2  linear  equations. 


C. 


(x,y) 


( k,k)=  Z  ajC(xy)(j,k) 


j=—M 

C  =  (  1-  I  ak)E[ffx)] 

k=-M 


(k  =  0,±1,  -±M) 


(15) 


and  E\fy(x)]  be  evaluated  as  E[f  (x)]  =  — - —  Z  fy  (x) . 

2M  + 1  k=-M 

C(x,y)(j,k )  is  adaptively  evaluated  using  the  neighboring 
pixels  around  (x,y)  in  [fy  (x)} , 
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<WM)  = 


1 


I  N(X,  y)  I 

(j,k  =  0  +!,••■  ±M) 


The  linear  MMSE  estimate  is  summarized  as  follows: 

1.  For  sample  (x,y),  evaluate  the  local  statistics  E\fy(x)] 
and  C(xyj(k,j)  using  the  estimates  {/* (x)}  in  the  local 


region.  Solve  (15)  to  get  optimal  { }  and  C.  Use  (6)  to 
compute  the  composite  estimate  of f(x,y). 

2.  Apply  the  same  procedure  to  every  sample  in  f. 

The  linear  MMSE  estimator  combines  the  estimate 
results  by  partial  information  in  an  optimal  way  such  that 
the  estimation  MSE  is  minimized.  This  estimator  imposes 
spatially  adaptive  filtering  on  signals  interpolated  with 
partial  information  and  subsequently  assures  the  intensity 
consistency  along  the  locally  enhanced  edges. 


4.  EXPERIMENT  RESULTS 


We  applied  the  linear  composite  MMSE  estimation 
algorithm  to  image  interpolation.  Compared  with  bilinear 
and  bicubic  schemes,  the  MMSE  estimation  scheme 
produces  clearer  boundaries  in  the  interpolated  images.  In 
addition  to  sharpening  strong  edges,  the  linear  estimation 
step  assures  the  intensity  consistency  along  edges  and 
avoids  unnecessary  artifacts  around  enhanced  edges. 

Figure  4  shows  the  example  of  a  subregion  from  the 
interpolated  Lena  image  by  different  methods,  where  the 
image  is  interpolated  to  four  times  the  original  size.  Since 
the  wavelet  coefficient  synthesis  deals  with  sparsely 
located  strong  edges,  as  we  expect,  the  improvement  of 
the  visual  quality  is  more  obvious  in  the  regions 
containing  strong  edges  than  in  the  regions  with  many 
weak  edges  such  as  in  textures.  In  Table  1,  we  show  an 
example  of  the  PSNR  improvement  by  the  proposed 
estimation  approach  with  M  set  to  1.  The  numbers  were 
obtained  from  interpolating  the  smoothed  Lena  image  to  2 
and  4  times  the  size.  The  full  size  (512x512)  image  was 
first  smoothed  by  the  lowpass  filter  and  downsampled. 
Then  the  smoothed  image  was  interpolated  to  full  size  and 
compared  with  the  original  image.  A  4~5dB  gain  is 
achieved  by  the  proposed  estimation  algorithm. 

5.  CONCLUSIONS 


In  this  work,  we  used  signal  estimation  techniques  to 
solve  the  interpolation  problem.  We  proposed  a  linear 
composite  MMSE  estimator  to  recover  the  fine -resolution 
image  from  the  low-resolution  samples.  An  edge  model 
based  nonlinear  estimator  synthesizes  the  important 
detailed  coefficients  with  partial  information  obtained 
from  low-resolution  samples.  A  standard  linear  estimator 
then  minimizes  the  MSE  of  the  recovered  image  using 
knowledge  of  local  statistics  of  previously  estimated 
signals.  Subsequently,  strong  edges  are  enhanced  by  the 


model  based  coefficient  synthesis  while  the  intensity 
consistency  along  edges  is  maintained  by  the  subsequent 
linear  estimator.  The  interpolation  results  demonstrate  the 
strength  of  the  statistical  approach  in  achieving  better 
image  quality  and  better  signal  estimation.  Further  study 
includes  extending  the  statistical  analysis  to  the  problem 
of  multiscale  image  modeling. 


(c) 


Figure  4.  Images  produced  by  different 
interpolation  methods,  (a)  Bilinear;  (b)  Bicubic; 
(c)  Linear  composite  MMSE  estimator. 


Method 

x  2  (dB) 

x  4  (dB) 

Bilinear 

28.84 

23.35 

Bicubic 

29.11 

23.26 

Proposed 

33.78 

28.56 

Table  1.  PSNR  comparison  of  different 
interpolation  methods 
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